function f=NUFFT2D(alpha,omega,precision,M)

%
% Yoel Shkolnisky, January 2010.

if nargin<3
    precision='single';
end

[b,m,q]=NUFFTParameters(precision);

if length(size(omega))~=2
    error('omega must be a 2 dimensional array')
end

if size(omega,2)~=2
    error('omega must be an array with 2 columns');
end
mM=m*M;
alpha=alpha(:);
n=numel(alpha);
mu=round(omega.*m);
if size(omega,1)~=n
    error('alpha and omega must have the same number of rows');
end

% All precomputations are borken into 1D ones.
Px=zeros(n,q+1,precision);
Py=zeros(n,q+1,precision);
for j=-q/2:q/2
    tmp1=(m.*omega(:,1)-(mu(:,1)+j)).^2;
    tmp2=((m.*omega(:,2)-(mu(:,2)+j)).^2);
    Px(:,j+q/2+1)=exp(-tmp1/(4*b))./(2*sqrt(b*pi));
    Py(:,j+q/2+1)=exp(-tmp2/(4*b))./(2*sqrt(b*pi));
end

if strcmp(precision,'single')
    tau=NUFFT2DauxiliaryMEXsingle(n,M,m,q,mu,Px.',Py.',alpha);
elseif strcmp(precision,'double')
    tau=NUFFT2DauxiliaryMEXdouble(n,M,m,q,mu,Px.',Py.',alpha);
else
    error('precision must be single or double');
end

clear Px Py

T=fftshift(ifftn(ifftshift(tau)));
T=T*mM*mM;
low_idx_M=-ceil((M-1)/2);
high_idx_M=floor((M-1)/2);
idx=low_idx_M:high_idx_M;
E=exp(b*(2*pi*idx./(mM)).^2);
offset1=ceil((mM-1)/2);
offset2=offset1+low_idx_M+1;
E2=reshape(E'*E,[M M]);
f=T(offset2:offset2+M-1,offset2:offset2+M-1).*E2;